/*
DenseMatrix.h
Written by Matthew Fisher

a dense matrix structure.
*/

template<class T>
class DenseMatrix
{
public:
    DenseMatrix();
    explicit DenseMatrix(UINT Size);
    explicit DenseMatrix(UINT RowCount, UINT ColCount);
    DenseMatrix(const DenseMatrix<T> &M);

    ~DenseMatrix();

    //
    // Memory
    //
    void Allocate(UINT RowCount, UINT ColCount);
    void Clear(T Value = 0.0);
    void FreeMemory();
    DenseMatrix<T>& operator = (const DenseMatrix<T> &M);

    //
    // Files
    //
    void WriteToStream(ostream &os, char delimiter) const;

    //
    // Accessors
    //
    __forceinline T* operator [] (int Row)
    {
        return &_data[Row * _cols];
    }
    __forceinline const T* operator [] (int Row) const
    {
        return &_data[Row * _cols];
    }
    __forceinline T& operator()(int Row, int Col)
    {
        return _data[Row * _cols + Col];
    }

    __forceinline T operator()(int Row, int Col) const
    {
        return _data[Row * _cols + Col];
    }
    __forceinline T& Cell(int Row, int Col)
    {
        return _data[Row * _cols + Col];
    }

    __forceinline T Cell(int Row, int Col) const
    {
        return _data[Row * _cols + Col];
    }

    __forceinline UINT RowCount() const
    {
        return _rows;
    }

    __forceinline UINT ColCount() const
    {
        return _cols;
    }

    //
    // Operators
    //
    void operator += (const DenseMatrix &M);
    
    //
    // In-place modifiers
    //
    void Identity();
    void Identity(UINT Size);
    void InvertInPlace();

    //
    // Query
    //
    bool Square() const;
    bool ElementsValid();
    T Determinant() const;
	T DeterminantLaplaceExpansion() const;
	void MinorMatrix(UINT Row, UINT Col, DenseMatrix<T>& result) const;
    void ExtractRow(UINT Row, vector<T> &result);
    void ExtractCol(UINT Col, vector<T> &result);
    vector<T> ExtractRow(UINT Row);
    vector<T> ExtractCol(UINT Col);
    
    //
    // Linear algebra
    //
    DenseMatrix<T> Inverse() const;
    DenseMatrix<T> Transpose() const;
    void LUDecomposition(DenseMatrix<T> &L, DenseMatrix<T> &U);
    void LUSolve(vector<T> &x, const vector<T> &b);
    bool EigenSystem(vector<T> &Eigenvalues, DenseMatrix<T> &Eigenvectors) const;
    bool EigenSystem(T *Eigenvalues, T **Eigenvectors) const;
    string EigenTest(const vector<T> &Eigenvalues, const DenseMatrix<T> &Eigenvectors) const;
    string EigenTest(const T *Eigenvalues, const T **Eigenvectors) const;

    bool EigenSystemVictor(vector<T> &Eigenvalues, DenseMatrix<T> &Eigenvectors) const;
    bool EigenSystemTNT(vector<T> &Eigenvalues, DenseMatrix<T> &Eigenvectors) const;
    
    //
    // Conversion
    //
    /*void ToColumnVector(vector<T> &result)
    {
        Assert(_cols == 1, "ToColumnVector called with more than 1 column");
        result.ReSize(_rows);
        for(UINT RowIndex = 0; RowIndex < _rows; RowIndex++)
        {
            result[RowIndex] = Cell(RowIndex, 0);
        }
    }*/

    //
    // Static helper functions for functional-style sparse matrix manipulation.  Since the copy constructor
    // for a dense matrix can be costly, this approach may be more efficient that using operator overloading.
    //
    static void Multiply(DenseMatrix<T> &result, const DenseMatrix<T> &left, const DenseMatrix<T> &right);
    static void MultiplyMMTranspose(DenseMatrix<T> &result, const DenseMatrix<T> &M);
    static void Multiply(vector<T> &result, const DenseMatrix<T> &left, const vector<T> &right);
    static void Multiply(T *result, const DenseMatrix<T>&left, const T *right);
    static void MultiplyInPlace(DenseMatrix<T> &result, T right);
    static DenseMatrix<T> OuterProduct(const vector<T> &A, const vector<T> &B);
    static T CompareMatrices(const DenseMatrix<T> &left, const DenseMatrix<T> &right);

    static vector<T> LinearRegression(const DenseMatrix<T> &X, const vector<T> &Y);

private:
    UINT _rows, _cols;
    T *_data;
};

//
// Overloaded operators
//
template<class T> ostream& operator << (ostream &os, const DenseMatrix<T> &m);

template<class T> vector<T> operator * (const DenseMatrix<T> &left, const vector<T> &right);
template<class T> DenseMatrix<T> operator * (const DenseMatrix<T> &left, const DenseMatrix<T> &right);
template<class T> DenseMatrix<T> operator * (const DenseMatrix<T> &left, T right);
template<class T> DenseMatrix<T> operator * (T left, const DenseMatrix<T> &right);
template<class T> DenseMatrix<T> operator + (const DenseMatrix<T> &left, const DenseMatrix<T> &right);
template<class T> DenseMatrix<T> operator - (const DenseMatrix<T> &left, const DenseMatrix<T> &right);

template<class T> OutputDataStream& operator << (OutputDataStream &S, const DenseMatrix<T> &M)
{
    UINT rowCount = M.RowCount();
    UINT colCount = M.ColCount();
    S << rowCount << colCount;
    UINT totalEntries = rowCount * colCount;
    const T *data = M[0];
    for(UINT entryIndex = 0; entryIndex < totalEntries; entryIndex++)
    {
        S << data[entryIndex];
    }
    return S;
}

template<class T> InputDataStream& operator >> (InputDataStream &S, DenseMatrix<T> &M)
{
    UINT rowCount, colCount;
    S >> rowCount >> colCount;
    M.Allocate(rowCount, colCount);
    UINT totalEntries = rowCount * colCount;
    T *data = M[0];
    for(UINT entryIndex = 0; entryIndex < totalEntries; entryIndex++)
    {
        S >> data[entryIndex];
    }
    return S;
}

#include "DenseMatrix.inl"